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SUMMARY 

A new efficient relaxation scheme In conjunction with a multigrid method 
Is developed for the Euler equations. The LU-SSOR scheme Is based on a cen- 
tral difference scheme and does not need flux splitting for Newton Iteration. 
Application to transonic flow shows that the new method surpasses the perform 
ance of the LU Implicit scheme. 


INTR ODUCTION 

Recently several Implicit schemes have been successfully developed In 
conjunction with a multigrid method for steady-state solution of the unsteady 
Euler equations (refs. 1 to 3). Although the alternating direction Implicit 
scheme could be Improved to achieve the expected efficiency of the multigrid 
method In two-dimensions (ref. 1), Its Inherent limitations In three- dimensions 
suggested alternative approaches (ref. 2). An alternative Implicit scheme 
which Is stable In any number of space dimensions was based on LU factorlza 
tlon. The LU Implicit scheme was proved to be robust and efficient for high 
Mach number flows as well as transonic flows (ref. 4). It was also shown that 
a symmetric Gauss-Seldel relaxation method for solving the unfactored Implicit 
scheme was a variation of the LU Implicit scheme. 

The Newton Iteration method has been a subject of Investigation for solu- 
tion of the steady Euler equations (refs. 5 to 7). Because of the rapid 
growth of the operation count with the number of mesh cells, the system was 
solved Indirectly. Jespersen (ref. 5) and Hemker and Spekreljse (ref. 6) used 
the symmetric Gauss-Seldel method while MacCormack (ref. 7) applied the line 
Gauss-Seldel method to the Navler-Stokes equations. In this paper an efficient 
relaxation scheme In conjunction with a multigrid method Is developed for 
approximate Newton Iteration. The new LU-SSOR scheme needs scalar diagonal 
Inversions while the Gauss-Seldel method or the LU Implicit scheme require 
block matrix Inversions. It Is desirable that the matrix should be diagonally 
dominant to assure the convergence of a relaxation method. The new method 
based on a central difference scheme achieves this without flux splitting which 
substantially Increases the computational work per cycle. 




GOV ER NING E QUATIONS 


The Euler equations are obtained from the Navler- Stokes equations by 
neglecting viscous terms. Let p, u, v, E, H, and p be the density, Car- 
tesian velocity components, total energy, total enthalpy, and pressure, and let 
x and y be Cartesian coordinates. Then for a two-dimensional flow these 
equations can be written as 


aw aF + 3G 
at + ax + ay 


0 


( 1 ) 


where W Is the vector of dependent variables, and F and G are convective 
flux vectors 


W = (p,pu,pv,pE) T 

2 T 

F = (pu.pu 4- p,pvu,u(pE 4- p)) (2) 

2 T 

G = ( pv , puv , pv 4- p,v(pE 4- p)) 

The pressure Is obtained from the equation of state 

P = p(y - 1){e - 2 ( u? + v? )} (3) 

These equations are to be solved for a steady state aW/at = 0 where t 
denotes time. 


S EMI - DISCRETE F INITE VOLU ME ME1 HOD 

A convenient way to assure a steady state solution Independent of the time 
step Is to separate the space and time discretization procedures. In semi- 
discrete finite volume method one begins by applying a semi-discretization 
in which only the spatial derivatives are approximated. The use of a finite 
volume method for space discretization allows one to handle arbitrary geom 
etrles and helps one to avoid problems with metric singularities that are 
usually associated with finite difference methods. The scheme reduces to a 
central difference scheme on a Cartesian grid, and Is second order accurate in 
space provided that the mesh Is smooth enough. It also has the property that 
uniform flow Is an exact solution of the difference equations. 

NONLINEAR ADAPTIVE DISSIPATION 

In typical calculation of flow with discontinuities by a central differ- 
ence scheme, wiggles appear In the neighborhood of shock waves where pressure 
gradient is severe. In order to suppress the tendency for spurious odd and 
even point oscillations, and to prevent unsightly overshoots near shock waves, 
the scheme Is augmented by artificial dissipative terms. The dissipative term, 
which is constructed so that It is of third order In smooth regions of the 
flow, Is explicitly added to the residual. For the density equation, for 
example, the dissipation has the form 


2 



d Hl/2,j ~ d i-l/2,j * d i,jH/2 ‘ d i,j-l/2 

where 


.( 2 ) 


d Hl/2,j = e i+l/2,j (p Hl,j • p 1,j ) 


.( 4 ) 


e i+l/2,j (p U2,j 3p Ul,j ' 3p i,j p iU : 


( 4 ) 


Let S be the cell area which is equivalent to the inverse of the determinant 
of transformation Jacobian. Both coefficients include a normalizing factor 
Sf^i /2 proportional to the length of the cell side, and j 

also made proportional to the normalized second difference of the pressure 


'•i 'l p Hl.j * 2p i,j * p 1-l.j 


( 5 ) 


in the adjacent cells. The third order terms provide background damping of 
high frequency modes. The first order terms are needed to control oscillations 
in the neighborhood of shock waves, and are turned on by sensing strong pres 
sure gradients in the flow. The dissipative terms for the other equations are 
constructed from similar formulas with the exception of the energy equation 
where the differences are of pH rather than pE. The purpose of this is to 
allow a steady state solution for which H remains constant. Increasing the 
amount of artificial viscosity improves the rate of convergence although too 
much dissipation can hurt it. However, it is desirable to make the amount be 
as small as possible in order not to degrade the accuracy of solution. Typical 
amount of the third order terms is almost negligible when compared to the 
physical viscosity. 


J.U-SSOR SCHEME 

A prototype implicit scheme for a system of nonlinear hyperbolic equations 
such as the Euler equations can be formulated as 


W n+1 = W n - 6 At(D F(W n+1 ) * 0 G(W n+1 )) (1-0) At( D F ( W n ) * U G(W n )) (6) 

X y x y 


where D x and D y are difference operators that approximate a/dx and 
a/ay. Here n denotes the time level. In this form the scheme is too expert 
sive, since it calls for the solution of coupled nonlinear equations at each 
time step. Let the Jacobian matrices be 


aF n aG 

" aw, ' aw 


(/) 


and let the correction be 


<sw - w nf | w n 
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The scheme can be linearized by setting 

F(W n+1 ) = F(W n ) + A6W + 0(IISWll 2 ) 

G(W"+1) = G( W n ) + B6W t 0(II6WII 2 ) 
and dropping terms of the second and higher order. 

This yields 

[ItB At( D X A f DyB) ]<5W t At R = 0 (8) 

where R Is the residual 

R = D x F(W n ) + D y G(W n ) 

If B = 1/2 the scheme remains second order accurate In time, while for other 
values of B the time accuracy drops to first order. 

The unfactored Implicit scheme equation (8) produces a large block banded 
matrix which Is very costly to Invert and requires huge storage. If B = 1 the 
scheme reduces to a Newton Iteration In the limit At -> °°. 

( D X A * 0 y B) 6W * R = 0 (9) 

A diagonally dominant form of equation (9) 


(DA + + DV 4- dV t- oV) 6W + R = 0 
x x y y 


( 10 ) 


can be written as 


A 1j 6W 1j " A 1-l,j 6W 1-U f A mj 6W H1,3 ' A ij 6W ij 


* B ij 6W 1j ~ B 1,j-l 6W U-l * B 1 ,j+l 4W i,j4l " B 1j AW 1j + R ij = 0 


(ID 


By simulating It with backward and forward relaxation sweeps, we obtain the 
symmetric successive over-relaxation (SSOR) method, which can be written in two 
steps as 


(A 1j - V 4W 1i f A 1*l J 4W U1J * < e ,3 - V 4W 1j * * "u 


4- R . . ^ 0 


followed by 


( 12 ) 


(A U - A i j )4W lJ ■ A lu“l-i.)' a h 1 j<,,j * Kj - 


l >R iJ i0 (13) 
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where D and D~ are backward difference operators and D* and D* are 

x y x y 

forward difference operators. Here, two-point operators are used for steady 

flow calculations. A + , A - , B + , and B~ are constructed so that the 

eigenvalues of "+" matrices are non-negative and those of matrices are 

non-positive. 

A+ = 2 (A + r A I} ’ A ~ = \ (A ' r A X) 

B f = \ (B t r B I), B' = \ (B - r B l) (14) 

where 

r A >. max (| \ A I), r g > max (| X g | ) (15) 

Here, and \q represent eigenvalues of Jacobian matrices. 

Subtract equation (12) from equation (13) to get 

(fl 1j ‘ A 1J> 4W 1J " *1-1 , j * <B 13 B 1j )4 “lj B 1J-1 4W i,J-l 

■ (A 1J " V 4l Sj * Kj - VKj (16) 

This may be written as 

( D~A + 4- D B + - A - B") 6W = (A + t B* - A' - B _ ) 6W* (17) 

x y 

where 

6W* = (DV + dV f A + t B*) ‘V R) (18) 

x y 

If we take and matrices as given In equation (14), then 

A + - A = r ft l, B f - B~ = r B l 

Thus equation (17) becomes the LU- SSOR scheme for approximate Newton iteration 
(D x A + + D~B + - A" - B') (dV + dV + A + 4- B*) <5W = - ( r r r g ) R (19) 

The equation (19) can be Inverted in two steps. 


MUL TIGRIO METHOD 

In order to adapt the LU-SSOR scheme for a multigrid algorithm, auxiliary 
meshes are Introduced by doubling the mesh spacing. Values of the flow vari- 
ables are transferred to a coarser grid by the rule 
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( 20 ) 


^ S h G) h 
" 2h = S 2h 

where the subscripts denote values of the mesh spacing parameter, S is the 
cell area, and the sum is over the four cells on the fine grid composing each 
cell on the coarser grid. The rule conserves mass, momentum, and energy. The 
solution on a coarse grid is updated as follows. 

(1) Calculate the correction and update the solution on the fine grid 

(2) Transfer the values of the variables to the coarse grid 

(3) Collect the residual on the fine grid for the coarse grid. A forcing 
function is then defined as 


P 2l>'K R 2h < 2 ’> 

where R is the residual. Superscripts c and t mean the collected and the 
transferred values respectively. The residual on the coarse grid is given by 

R 2h = R 2h ‘ P 2h ' K < 22 > 

(4) Calculate the correction and update the solution on the coarse grid. 
For the next coarser grid the residual is recalculated as 


R 4h * R 5h * P 4h - * < R 2h * P 2h> 
Similarly, the residual for the next grid is 


* < R 2h - K 


R 2h> 


R Sh - * £ < R 2 U h * £ R 


R 2h* 


<u 


(23) 


(24) 


where the superscript u means the updated value. On the first coarse grid, 

is replaced by £R^ with the result that the evolution on the coarse 

grid is driven by the residuals on the fine grid. The evolution on the next 
coarser grid is driven by an estimate of what the fine grid residuals would 
have been as a result of the correction on the first coarse grid. The process 
is repeated on successively coarser grids. 


(5) Finally, the correction calculated on each grid is passed back to the 
next finer grid by bilinear interpolation. 

Since the evolution on a coarse grid is driven by residuals collected from 
the next finer grid, the final solution on the fine grid is independent of the 
choice of boundary conditions on the coarse grids. The surface boundary con- 
dition is treated in the same way on every grid, by using the normal pressure 
gradient to extrapolate the surface pressure from the pressure in the cells 
adjacent to the wall. Values are extrapolated to the fictitious cells Inside 
the body surface for the second difference dissipation on the coarse grids. 

The far field conditions can either be transferred from the fine grid, or 
recalculated. 
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RESULTS 


The new LU-SSOR scheme In conjunction with a multigrid method Is applied 
to transonic airfoil calculation. The test case Is the NACA 0012 airfoil at 
Mach 0.8 and 1.25° angle of attack. Figures 1 and 2 show the plot of Mach 
number contours and the surface pressure distribution respectively. These 
solutions are obtained on a 128 by 32 C-mesh. The convergence histories are 
shown In figures 4 and 6. A four-level multigrid Is used without grid 
sequencing. The lift history Is shown In figure 4 while the density residual 
Is shown In figure 6. The results obtained by the LU Implicit scheme are pre- 
sented In figures 3 and 5 for comparison. As the results show, the convergence 
rate of the LU-SSOR scheme Is about 30 percent faster than that of the LU 
Implicit scheme. Moreover, the computational work for the LU SSOR scheme Is 
about 30 percent less than that for the LU Implicit scheme. The overall com- 
putational work Is reduced by a factor of two. The convergence for the engi- 
neering accuracy Is usually achieved In less than 10 CPU seconds with the Cray 
XMP computer. 
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Figure 6. - Convergence history (present method). 
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